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We present an efficient integral equation approach to solve the heat equa- 
tion, Wt(x) — Am(x) = -F(x, t), in a two-dimensional, multiply connected 
domain, and with Dirichlet boundary conditions. Instead of using integral 
equations based on the heat kernel, we take the approach of discretizing in 
time, first. This leads to a non-homogeneous modified Helmholtz equation 
that is solved at each time step. The solution to this equation is formulated 
as a volume potential plus a double layer potential. The volume potential 
is evaluated using a fast multipole-accelerated solver. The boundary condi- 



lO \ tions are then satisfied by solving an integral equation for the homogeneous 

modified Helmholtz equation. The integral equation solver is also acceler- 
ated by the fast multipole method (FMM). For a total of N points in the 

£->! . discretization of the boundary and the domain, the total computational cost 

per time step is O(N) or 0(N log N). 

Keywords: Fast multipole method, Gaussian quadrature, Modified 
Helmholtz equation, integral equations, the heat equation. 
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1. Introduction 

Integral equation methods offer an attractive alternative to conventional 
finite difference and finite elements for solving partial differential equations 
that arise in science and engineering. They offer several advantages: com- 
plex physical boundaries are easy to incorporate, the ill-conditioning associ- 
ated with directly discretizing the governing equation is avoided, high-order 
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accuracy is easier to attain, and far-field boundary conditions are handled 
naturally. Nevertheless, integral equation methods have been slow to be 
adopted for large-scale problems, and the reason for this is clear. For non- 
homogeneous problems, unknowns are distributed at all points of the domain. 
On an N x N mesh, integral equation formulations lead to dense N 2 x iV 2 
matrices. This situation is further exacerbated if the problem is time depen- 
dent. 

In recent years, fast algorithms have been developed for integral equa- 
tions in a variety of settings. An early example is the work by Greenbaum et 
al. in [6j which formulates integral equation methods for Laplace's equation 
in multiply connected domains. The solution to these integral equations was 
accelerated by the original version Greengard and Rokhlin's fast multipole 
method (FMM) [9|, |lOj. For iV points in the discretization of the domain, 
the solution to the integral equation requires only O(N) operations. Subse- 
quently, similar tools have been developed for other elliptic partial differential 
equations, such as Poisson's equation [5j, the Stokes equations [8], and the 



modified Helmholtz equation^, [13|. Nishimura provides a comprehensive 



review of fast multipole accelerated integral equation methods in [16 . 

For time-dependent problems, there is added complexity. There is recent 
work by Li and Greengard in [14], [l5| on developing fast solvers for the heat 
equation. Their methods use an integral equation formulation based on the 
heat kernel, and while robust, are highly complicated. To date, the full set 
of tools required for solving the heat equation in interior or exterior domains 
with complex boundaries has not been constructed. 

We present an alternative approach to solving the heat equation and other 
time-dependent partial differential equations that rely on available fast algo- 
rithms for elliptic problems. The temporal discretization of the heat equation 
gives rise to the modified Helmholtz, or linearized Poisson-Boltzmann, equa- 
tion 

M(x)-tt 2 A«(x) = 5(x,i), (1) 

where a 2 = O(At). This equation seems to have received little attention 
in the literature, particularly in this context. Hsiao and Maccamy formulate 



integral equations of the first kind in |12[ for this equation. In [4], Cheng et al 



present a fast direct solver for (0Q) in two dimensions on the unit square. The 
solution is expressed as a volume potential, and the direct solver is accelerated 
using a new version of the fast multipole method [7|, |ll| . The solver is fully 
adaptive and the computational costs are comparable to those of FFT based 




Figure 1: The physical domain f2 (in gray) is embedded in the unit square D. The physical 
outer boundary is denoted by Tq and the interior boundary is T±. 



methods. In [13], we present fast, well-conditioned integral equation methods 



to solve the homogeneous equation in unbounded and bounded multiply- 
connect domains, with Dirichlet or Neumann boundary conditions. What we 
present here is a preliminary study on coupling the methods discussed in J4J 
to those discussed in [l3| in order to solve practical problems in complicated 
domains. 

The paper is organized as follows: In Section 2, we discuss the temporal 
discretization of the heat equation and how this gives rise to the modified 
Helmholtz equation. We outline the corresponding potential theory; decou- 
ple the solution into a volume potential and double layer potential with an 
unknown density, and then formulate an integral equation for the Dirichlet 
problem in a bounded multiply connected domain. In Section 3, numerical 
methods for the evaluation of the volume potential and the integral equation 
are discussed. Numerical examples are presented in Section 5. 



2. Potential theory and the integral equation formulation 

Consider the isotropic inhomogeneous heat equation in a bounded domain 
ft in R 2 : 

ut(x,t) - Au(x,t) = F(x,t), £>0, (2) 

with Dirichlet boundary conditions: 



u(x, t) = /(x, t), x e r, t > 0, 



and initial conditions 

w(x, 0) = w (x), xGfi. 

One approach to solving this boundary value problem would be base an 
integral equation approach on the fundamental solution for the heat equation. 



Fast algorithms are developed for just such an approach in |14j . |15fl. We 
approach from a different starting point by discretizing (T5]) in time, first, and 
then reformulate in terms of an integral equation. 

In order to prevent a severe time-step restriction, linearly implicit or 
IMEX schemes [3J are generally used for marching in time. In such schemes, 
the diffusive term is treated implicitly, while the remainder is treated explic- 
itly. Regardless of the details of the particular choice of IMEX scheme, the 
temporal discretization of (J2J) yields 

u N+1 - a 2 Au N+1 = B{u N ,u N -\---), xell, 

^ +1 = /(x,t), xer, (3) 

u° = ito(x), at t = 0, 

where t = (N + l)At is the current time. The simplest such scheme is the 
first-order backward Euler method, which is of the form ([3]) with 

a 2 = At, B = u N + AtF N . 

A second-order method is the extrapolated Gear method, for which 

a 2 = -At, B = -u N - V" 1 + -AtF N - -AtF N ~\ 
3 ' 3 3 3 3 

In order to solve ([3]) by means of integral equations, we first represent the 
solution at each time step as 

[/(x) = £/ p (x) + £/ /l (x), 
U p - a 2 AU p = B, x G Q, (4) 

U h - a 2 AU h = 0, xell, 

u h = 0( X ), xer, ! '' 



where 

and U h satisfies 



where g(x) = /(x) — t/ p (x), x G T. 
The solution to (J3J is represented as 

C/*(x) = f B{y)G{y-*)dA y , (6) 

where the free space Green's function G(x) for the operator 1 — a 2 A is given 
by the zeroth-order modified Bessel function of the second kind, 

We seek the solution to (jSJ) in the form of a double layer potential 

where cr(y) is the value of an unknown density at the boundary point y, 
and d/dvy represents the outward normal derivative at the point y. The 
unknown density a is found by solving an integral equation derived to ensure 
the boundary conditions in (jSJ) are satisfied. 



In [13j, we derive integral equation formulations to solve (j5|) in bounded 
and unbounded multiply-connected domains. Here, we summarize the results 
for the bounded case. Observing that K (z) ~ —log(z) as z — > [1], we can 
determine the jump relations of the double layer potential as these are well 
known for the logarithmic kernel. Thus for any point x on the boundary T, 

x' € SI 

Substituting the above into the boundary condition in (J5J) results in the 
integral equation 

' *(*) + 77^ / 4-Ko ( lZ ^-) <r(y) ds(y) = g(x). (7) 



2a 2 2na 2 J T du y \ a 

The kernel in the above is continuous along T, 

lim jr—K I I = --re(x), 

x, yer 

where k(x) denotes the curvature of T at the point x, (here, we have again 
used the logarithmic behaviour of K (z) as z — > 0). In addition, ([7|) has 
no nontrivial homogeneous solutions. Therefore, by the Fredholm alterna- 
tive, (JZ|) has a unique solution for any integrable data /(x). 



3. Numerical methods 

We now discuss the numerical methods for solving the modified Helmholtz 
equation ([3]). In Section 3.1, we first outline the methods involved in evalu- 
ating the volume potential flS]). In Section 3.2, we next discuss the coupling 
of the volume potential to the double layer potential through the boundary 
conditions and solving the corresponding integral equation (J7J). The work 



presented in this section is discussed in more detail in [131 ]. While the FMM 
plays an integral role in our numerical methods, its implementation is stan- 
dard and has been discussed at length in other works. For details on the 
FMM, we refer the reader to the work in [ly, |9] for the original FMM, 11 
for the new FMM algorithm based on exponential expansions to represent the 
far-field interactions, and |4| for the FMM applied to the modified Helmholtz 
equation in two dimensions. 

3.1. The volume potential 

In |4j, methods for the rapid evaluation of the volume potential for the 
modified Helmholtz equation on the unit square were developed. These meth- 
ods use adaptive mesh refinement, and the evaluation is accelerated using the 
new version of the FMM. The total computational cost is comparable to that 
of an FFT-based fast solver. In order to employ these methods, we must first 
extend the right-hand side of (|2j), which is defined in Q, to be defined through- 
out D. One naive way to define an extension -B(x, t) of B(x,t) to D is the 
following: 

B(x,t), xGfi 



B(x,t) "\ o, xefl\fi. 

However, the discontinuity across T would likely cause excessive grid re- 
finement in the vicinity of the boundary, and consequently slow down the 
evaluation (c.f. example 4.4 in |5J; up to O(10 5 ) grid points are needed to 
solve Au = 1 inside a disk and Am = outside). For multiple-component 
boundaries, this would simply become too expensive. For the purposes of 
this present work, we consider problems with simple boundary conditions, 
only, so that B can be extended continuously from Q to D. (We save the 
consideration of more general-purpose methods to extend B for future work.) 
For example, if we are solving the homogeneous heat equation (F(x., t) = 0, 
then -B(x) consists of a linear combination of the solution w(x, t) at previous 
time steps. If the boundary conditions are constant, u = Co on To, u = C\ 



on I\, say, then we can define a continuous -B(x, t): 



B(x,t), 


x e 0, 


Cq, 


exterior to r 


Ci, 


exterior to I\ 



We now use the methods presented in [4J to evaluate 

U p (x)= [ B(x,t)G(x-y)dA y , (8) 

J D 

and let £/ p (x) = £>(x) for xefl. 

The FMM uses an adaptive quad-tree structure in order to superimpose 
a hierarchy of refinement on the computational domain. The unit square D 
is considered to be grid level 0. Grid level / + 1 is obtained recursively by 
subdividing each square (or node) s at level / into four equal parts; these are 
called the "children" of s. Adaptivity is achieved by allowing different levels 
of refinement throughout the tree. The square s is refined if the error of a 
fourth order polynomial interpolating B is larger than some preset tolerance. 
We denote the childless nodes in the quad-tree as -Dj, i — 1, ■ ■ ■ , M, where 
M is the total number of such nodes. We assume we are given Bona cell- 
centered 4x4 grid for each Di. Thus, Nq = 16 x M is the total number of 
grid points in D. To obtain fourth-order accuracy, these 16 points are used 
to construct a fourth-order polynomial to B of the form 

10 
3=1 

where x* is the centre of D^ and {pi\ are the standard basis functions for 
polynomials up to order three (see [4j, |5|] for details on the approximating 
polynomial). Therefore, the evaluation of © is approximated by 



M r 10 

£>(x)^W G(x-y)^c}^(x-x' ; )^ 



=1 JD * 3=1 

Direct evaluation of this potential at all Np grid points would require O(N^) 
work. This expense is reduced to O(Nd) work, by using the FMM. 

Once we have obtained the values of U p at all Np points, we must next 
evaluate the volume potential on V in order to impose the correct boundary 



conditions in 
on each D*. 



). Again, we construct the approximating polynomial to U p 



10 



f/ p (x)^^q Pi (x-x'), xeA- 



Every point y G T belongs to a particular childless node Di. Determin- 
ing the value of U p at y is simply a matter of evaluating the appropriate 
approximating polynomial at that point. 

3.2. The integral equation 

We now discuss the numerical methods to solve ([7]). We assume each 
component curve Tk, k — 0, 1 is parametrized by y k (a), where a G [0,2-7r). 
Similarly, a k (a) refers to the restriction of the density a on IV We are given 
N points equispaced with respect to a. Thus the mesh spacing is h = 2n/N, 
and the total number of descretization points is iVr = 2N. Associated with 
each such point, denoted by y k , is an unknown density a k . 

In order to approximate the integral operator in j[7|) , we use hybrid gauss- 
trapezoidal quadrature rules developed by Alpert [2j which are tailored for 
integrands with logarithmic singularities. These quadratures are of order 
h p \og h. The order p determines the nodes u n and weights v n , n = 1, • • • , I, 
which are used for the quadrature within the interval a G [a, — ha, ctj + ha], 
on Tfc (a is also determined by p). Outside of this interval, the quadrature is 
essentially the trapezoid rule. Applying this quadrature to (J7J) yields 



^ h / 



1 N 



N+j-a 



m = 
m 7^ k 



n=l 



n=j+a 



1 ' 



• 2 -'yS) 



n = —l 



(9) 



where 



a \ a 1 y — x 



n y . 



The outward pointing normal n y at each point is obtained pseudospectrally. 
In the middle sum, we invoke periodicity of all functions on T^, or equiva- 
lently, j + N = j. In the final sum, we are required to know a at intermediate 
values to the nodal values. These are found through Fourier interpolation. 

Equation (J5J) is a linear system that is solved iteratively using the gener- 
alized minimum residual method GMRES [17] . The bulk of the work at each 
iteration lies in evaluating OH]) at the current solution update. If this was 
done directly, it would require O(N^) work. This evaluation can be reduced 
to O(Nt), again using the FMM (c.f. [13] for more discussion on the details 
of implementation). Since the number of iterations needed to solve a Fred- 
holm equation of the second kind to a fixed precision is bounded independent 
of N r , we can estimate the total cost of solving (J2J) by 

I(e)C{e)N r , 

where 1(e) is the number of GMRES iterations needed to reduce the residual 
error to e, and C(e) is the constant of proportionality in the FMM. 

4. Numerical results 

The algorithms described above have been implemented in Fortran. Here, 
we illustrate their performance on a variety of examples. The tolerance for 
the residual error in GMRES is set to 10~ 10 . All timings cited are for a Mac 
Pro 2.1 with two 3GHz Quad-Core Intel Xeon processors. 

EXAMPLE 1. In this example, we solve the forced heat equation on a 
domain for which an analytical solution is known. The domain Q is bounded 
between r , a circle of radius 0.4, and r 1; a circle of radius 0.1. Both T and 
Ti are centered at the origin. The forcing term is 

, ,n sin(20|x|) 

F(x,t) = 400cos(20|x|) + 20 — V ' u 

The following is the exact solution to the heat equation: 

u(x, t) = e~ xH [Fo(O.lA) J (A|x|) - J (0.1A)F (A|x|)] + cos(20|x|), 

where Jo is the Bessel function of the first kind of order zero, and Yq is 
the Bessel function of the second kind of order zero. We have chosen A so 



that the time-dependant term vanishes on both boundaries A 
boundary conditions, then, are 



10.244. The 



/(*) 



cos(8), 
cos(2), 



xG r , 
xe r\. 



We calculate solutions to the heat equation using the IMEX Euler and SBDF 
methods, up to t = l.Oe — 2. The results are summarized in Table [TJ 



At 
2.0 x 10~ 3 
1.0 x 10~ 3 
5.0 x 10~ 4 
2.5 x 10~ 4 



Error i 



Error2 
1.65 x 10" 3 
4.74 x 10" 4 
1.23 x 10" 4 
3.25 x 10~ 5 



1.37 x 10" 2 
7.09 x 10~ 3 
3.59 x 10~ 3 
1.78 x 10~ 4 



Table 1: Temporal Error using IMEX Euler and extrapolated Gear Method. Error i indi- 
cates the error for the IMEX Euler method, Error2 is for extrapolated Gear. 




Figure 2: The quadtree structure for Example 2. 

EXAMPLE 2. In this example, we solve the homogeneous heat equation 
in an elliptical region containing an off-center elliptical hole. The boundary 
conditions are 

'o, x g r , 
l, xeiY 



/(x) 



The quadtree structure for the solution homogeneous equation is shown in 
figure [2] The results are shown in figure [3j 
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Time =0.01 



fe D: ^ D: ^ 



Figure 3: The solution to Example 2. N r = 1024, N D = 65536. Total CPU time for 100 
time steps about 30 minutes. 

EXAMPLE 3. In this example, we demonstrate that our methods can 
be applied to much more complex equations. Here, we solve the Allen-Cahn 
equation 



u t — eAu = u(l - 


-A 


xe(], t > o, 


u(x, t) = 0, 




xer 


w(x, t) = 0, 




xeTi, 



where e = 10 -5 . We initialize the solution with random values uniformly 
distributed on [— |, |]. The general behaviour of solutions to the Allen-Cahn 
equation is well known: the stable stationary solutions are u = 1 and u = — 1 
and the solution exhibits coarsening towards these values. The presence of 
physical boundaries can create more complex patterns, as seen in figure |U 

5. Conclusions 

We have presented an investigation on coupling available fast algorithms 
for integral equations in order to investigate problems of scientific interest. 
We have shown that for time-dependent problems, doing a temporal dis- 
cretization first, and then reformulating the resulting elliptic equation as an 
integral equation, appears to have significant potential for handling compli- 
cated time dependent problems on complex domains. This paper is meant to 
be a preliminary investigation only; more work still has to be done to make 
this approach amenable to large-scale problems with a variety of boundary 
conditions. Specifically, we need to develop a more robust and general pur- 
pose method of extending functions defined on a complex domain onto a 
regular domain in a way that does not result in excess grid refinement in 
the vicinity of physical boundaries. Future applications will likely include 
developing integral equation methods for the incompressible Navier-Stokes 
and Euler equations. 
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Figure 4: The solution to Example 3, 100 time steps of size one are taken. The first 
picture shows the initial conditions, and the final picture is at t = 100. 
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